home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / jpl_c.zip / SIMPSON.C < prev    next >
Text File  |  1986-05-18  |  2KB  |  49 lines
  1. /* 1.0  03-05-85 */
  2. /************************************************************************
  3.  *                      Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, California 91009    1985        *
  6.  ************************************************************************/
  7.  
  8. #include "defs.h"
  9. #include "stdtyp.h"
  10. #include "errno.h"
  11.  
  12. /*----------------------------------------------------------------------*
  13.     This algorithm utilizes Simpson's extended rule, as published in
  14. HANDBOOK OF MATHEMATICAL FUNCTIONS, U. S. Department of Commerce, 
  15. National Bureau of Standards, Applied Mathematics Series 55, June 1964,
  16. formula 25.4.6 on page 886.  The error is (n*h**5/90)*f''''(x) for some
  17. x in (lowlim, uplim), where h = (uplim - lowlim)/2n.            */
  18.  
  19. /************************************************************************/
  20.     double
  21. simpson(func, lowlim, uplim, n)    /* return the integral of func from
  22.                    lowlim to uplim using Simpson's rule
  23.                    on 2n + 1 points.            */
  24. /*----------------------------------------------------------------------*/
  25.  
  26. double (*func)(), lowlim, uplim;
  27. {
  28.     double x, y, odds, evens, h, h2, ends;
  29.     int i;
  30.  
  31.     if (n <= 0)
  32.     {    errno = EDOM;
  33.         return 0.0;
  34.     }
  35.     odds = evens = 0.0;
  36.     ends = (*func)(lowlim) + (*func)(uplim);
  37.     h2 = (uplim - lowlim) / n;
  38.     h = h2 / 2.0;
  39.     x = lowlim + h;        /* sum upward on odd points, downward    */
  40.     y = uplim - h2;        /* on the even points.            */
  41.     for (i = 1; i < n; i++, x += h2, y -= h2)
  42.     {    odds += (*func)(x);
  43.         evens += (*func)(y);
  44.     }
  45.     odds += (*func)(x);
  46.     return (h * (ends + 4.0 * odds + 2.0 * evens) / 3);
  47. }
  48.  
  49.